Example data analysis workflow

Make up some data

import altair as alt
from IPython.core.display import Image, display
from bayes_window import generate_fake_spikes, fake_spikes_explore, BayesRegression
df, df_monster, index_cols, firing_rates = generate_fake_spikes(n_trials=20,
                                                                n_neurons=6,
                                                                n_mice=3,
                                                                dur=5,
                                                                mouse_response_slope=40,
                                                                overall_stim_response_strength=5)

Problem statement

There are three subjects with five neurons each. Some invervention is applied (pharmaceutical, stimulation, etc) that affects how often neurons are activated.

  • Neuron 4 of mouse 0 has the least effect

  • Neuron 0 mouse 2 has the most effect

  • Slow neurons are affected the most across all mice

We’ll consider interspike interval (ISI) as the measure of activity

Exploratory plots

# Make:
charts = fake_spikes_explore(df=df, df_monster=df_monster, index_cols=index_cols)

# Render:
charts[-3].properties(title='Change in firing rate (Hoverable, try it!)').display()
[chart.display() for chart in charts[:-1]];

Fully-Bayesian generalized linear model

Uses my homegrown bayes-window package for rapid prototyping

bw = BayesRegression(df=df,
                     y='isi',
                     treatment='stim',
                     condition=['neuron', 'mouse'], 
                     group='mouse',
                     detail='i_trial',
                     ).fit(dist_y='student',
                           num_warmup=1000,
                           n_draws=1000, 
                           num_chains=4,
                           add_group_intercept=True,
                           add_group_slope=True)

The effect on each neuron from each mouse is reconstructed approximately correctly. Faint band below is 95% highest density probability interval (HDPI), bright band 75% HDPI. Line overlays the mean for comparison

bw.plot(x=alt.X('neuron:N'), error_type='band').facet(column='mouse').properties(title='Posteriors. Also hoverable')

The overall firing rate for each mouse is also correctly reconstructed. Compare ISI means (lines connecting stim=0 and stim=1) to model intercepts (ticks):

bw.plot_intercepts().properties(title='Ticks=posteriors for intercepts, lines=raw ISI means ')

The overall effect for each mouse is correctly reconstructed:

slopes = bw.trace.posterior['slope_per_group'].mean(['chain', 'draw']).to_dataframe().reset_index()
chart_slopes = alt.Chart(slopes).mark_bar().encode(
    x=alt.X('mouse_:O', title='Mouse'),
    y=alt.Y('slope_per_group', title='Slope')
)
chart_slopes.properties(title='Quick and dirty estimate for group-level sloopes')

Out-of-sample prediction

In this synthetic dataset (simplified version for speed), gamma-distributed multilevel GLM (bw_gamma) works better for out-of-sample prediction than normally-distributed multilevel GLM (bw_normal) or frequentist mixed linear model (mlm):

display(Image('roc.png'))
display(Image('auc.png'))
../_images/neurons_example_16_0.png ../_images/neurons_example_16_1.png